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ABSTRACT 

We examine the possible effects of UV background fluctuations on the pixel flux power 
spectrum and auto-correlation function of the Lya forest due to a finite number of 
sources in an attenuating medium. We consider scenarios in which QSO sources dom- 
inate the contribution to the UV background. To estimate their contribution, we use 
the QSO luminosity functions from the 2dF and the Sloan Digital Sky Survey. We es- 
timate self-consistent values for the attenuation length at the Lyman edge through the 
Intergalactic Medium using Particle Mesh simulations of the Lya forest, normalised by 
the measured mean Lya flux. It is necessary to add the contribution from Lyman Limit 
Systems based on their measured statistical properties since the simulations are unable 
to reproduce the measured abundances of these systems, suggesting their formation 
may involve more complex hydrodynamical processes than simple collapse into dark 
matter halos. We examine the convergence properties of the flux power spectrum and 
auto-correlation function under differing assumptions regarding pressure-smoothing. 
Convergence to better than 10% in the flux power spectrum at proper wavenumbers 
k < 0.015 (kms~^)~^ requires a comoving resolution of 30 — 60 h^^ kpc in a comoving 
box of size 25 Mpc at z = 5, with the requirements becoming more severe at lower 
redshifts. The flux auto-correlation function does not converge to better than 10% on 
scales exceeding 3% of the box size. We find much better covergence properties for 
the relative effects of the UV background fluctuations on the flux power spectrum. 
The background fluctuations increase the large-scale power and suppress the power 
at intermediate scales. At z < 4, the effect is only at the few percent level. By z > 5, 
however, the large-scale power is boosted by more than 50%, suggesting the flux power 
spectrum may serve as a useful tool for distinguishing QSO-dominated UV background 
scenarios from those in which more abundant sources like galaxies dominate. 
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1 INTRODUCTION 

Due in large part to insights gained from hydrodynami- 
cal simulations of structure formation in the universe (Cen 
et al. 1994; Zhang, Anninos & Norman 1995; Hernquist 
at al 1996; Zhang et al. 1997; Bond & Wadsley 1997; The- 
uns, Leonard & Efstathiou 1998) the Lya forest, as mea- 
sured in high redshift Quasi-Stellar Object (QSO) spectra, 
has emerged as one of our premier probes of the universe 
at high redshift. The results of simulations assuming early 
homogeneous H i and He ii reionization reproduce the cu- 
mulative flux distributions and H i column density distri- 
butions measured from high resolution spectra to an accu- 
racy of a few percent (Meiksin, Bryan & Machacek 2001). 
More problematic are the predicted widths of the absorp- 
tion features, which appear to require additional sources of 



broadening, perhaps late He ii reionization, to match the 
measured widths (Theuns et al. 1999; Bryan & Machacek 
2000; Meiksin et al. 2001). 

Parallel to the development of the full hydrodynam- 
ical simulations of the Intergalactic Medium (IGM) are 
pseudo-hydrodynamical schemes using pure gravity (Petit- 
jean, Miicket & Kates 1995; Croft et al. 1998; Gnedin & Hui 
1998; Meiksin & White 2001). These simpler treatments are 
based on two key results from the full hydrodynamical com- 
putations. As a consequence of photoionization equilibrium: 
1. the gas follows a tight temperature-density relation for 
the absorbing material with the neutral hydrogen tempera- 
ture proportional to the baryon density to the power of ~ 0.5 
(Zhang et al. 1998; Schaye et al. 1999), and 2. the neutral 
hydrogen density closely traces the total matter density on 
the scales relevant to the forest (0.1 — 10/i~^Mpc) (Zhang 



© 0000 RAS 



2 A. Meiksin and M. White 



et al. 1998). The structure in QSO absorption thus tracks, 
in a calculable way, slight fluctuations in the matter density 
of the universe back along the line of sight to the QSO, with 
most of the absorption by the Lya forest arising from gas 
with ovcrdcnsitios of a few times the mean density. Although 
these techniques are not sufficiently accurate to be used for 
detailed modelling of the flux and absorption line proper- 
ties of the IGM (Meiksin & White 2001), they provide a 
useful means of investigating convergence questions and the 
possible role of additional physical effects on the Lya forest 
without the substantial overhead of the full hydrodynamical 
treatment. In particular, in recent years they have been in- 
creasingly relied on for predicting the flux power spectrum 
of the Lya forest (Meiksin & White 2001; Zaldarriaga, Hui 
& Tegmark 2001; Croft et al. 2002b; Seljak, McDonald & 
Makarov 2003). 

While much work has been done on the physics of the 
Lya forest and the cosmology which may be extracted, the 
role of fluctuations in the UV background on the absorp- 
tion properties of the Lya forest are largely unexplored 
(see Zuo 1992a,b; Fardal & ShuU 1993; Croft et al. 1999, 
2002a; Gnedin & Hamilton 2002). The effects of fluctua- 
tions in the UV background on the statistical properties of 
the Lya forest were considered previously by Zuo (1992a,b) 
and Zuo & Phinney (1993). Zuo & Bond (1994) introduced 
the use of the flux transmission auto-correlation function 
for measuring the clustering of the absorbing gas, although 
they did not consider the effect of fluctuations in the UV 
background apart from a local proximity effect. As we shall 
show, the effects of UV background correlations are most 
conspicuous only at high redshifts {z > 5), where the pre- 
cise level of absorption by the Lya forest is most difficult to 
measure because of the extremely low flux values, requiring 
very high signal-to-noise ratio spectra. 

Fardal & Shull (1993) explored the effects of UV back- 
ground fluctuations on the Lya forest along the linos of Zuo, 
but using Monte Carlo realisations of randomly distributed 
sources and clouds to estimate the effect of the background 
fluctuations. Croft ct al. (1999) estimated the flux power 
spectrum resulting from UV background fluctuations gen- 
erated by randomly distributed QSO sources in a uniform 
IGM, and found this to be small, at the few percent level, 
compared with the flux power spectrum from the Lya forest 
at 2 = 2.5. They found that allowing for QSO clustering had 
little additional effect. Croft et al. (2002a) accounted for the 
UV background fluctuations at z — 3 assuming the IGM was 
ionized by the gala^xies in their simulation, and concluded 
again the effect on the flux power spectrum is small. Finally 
Gnedin & Hamilton (2002) examined the effect of a fluctu- 
ating UV background at z = 4 using hydro simulations with 
pseudo-radiative transfer, assuming the reionization is ac- 
complished by galaxies in their simulation. They concluded 
the fluctuations of the UV background affect the flux power 
spectrum at the percent level, but cause an underestimate of 
the photo-ionization rate by 0(20%). A comparable boost 
in the photo-ionization rate is required if QSOs are the dom- 
inant source of the UV background at 2 > 5 (Paper I). 

In this paper, we investigate the effect of the fluctua- 
tions on the 2-pt flux distribution (the power spectrum and 
auto-correlation function) of the Lya forest. This extends 
our earlier work on the 1-pt distribution function reported 
in Meiksin & White (2003; hereafter Paper I). In the next 
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0.30 


0.70 
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0.70 


1.05 


0.97 


1.59 
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0.35 


0.65 


0.020 
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0.95 


0.88 


1.35 
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0.40 


0.60 


0.020 


0.55 


1.10 


0.90 


1.57 
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0.30 


0.70 


0.022 


0.70 


0.95 


0.92 


1.32 
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0.30 


0.70 


0.018 


0.67 


1.00 


0.90 


1.28 



Table 1. Parameters for the cosmological models. Hm is the total 

mass density parameter, Qy the vacuum energy density parame- 
ter, ^5 the baryonic mass fraction, h = i/o/100 km s^-' Mpc^-'^ , 
where Hq is the Hubble constant at z = 0, n tlie slope of tlie pri- 
mordial density perturbation power spectrum, erg the fluctuation 
normalization at 2 = in a sphere of radius &h~^ Mpc, and aj 
the fluctuation normalization at z = 3 flltered on the Jeans scale 
for gas at T = 2 X 10* K. 



section we describe the simulations used. In Section 3, we 
discuss our model for the UV background fluctuations. The 
results are presented in Section 4, and we summarise our 
conclusions in Section 5. In Appendix A wo discuss the nu- 
merical convergence tests we've performed, and we describe 
our estimates of the mean Lya flux in Appendix B. 



2 NUMERICAL SIMULATIONS 

To investigate the effects of UV background fluctuations on 
the power spectrum of the Lya forest, we have performed 
simulations of 5 different cosmological models within the 
ACDM family. In each case we used pure Particle Mesh 
(PM) dark matter simulations, mimicking the temperature 
of the gas using a polytropic equation of state and assum- 
ing gas and dark matter have the same spatial distribution. 
This has been shown (Petitjean, Miicket & Kates 1995; Croft 
et al. 1998; Gnedin & Hui 1998; Meiksin & White 2001), to 
produce results comparable to the full hydrodynamical sirrm- 
lations at the level of 10 per cent. The agreement is expected 
to be even better on the large scales of interest to this work. 
The parameters for the simulations are provided in Table 1. 
The models are all for a flat universe and are consistent 
with Lya forest constraints on matter fluctuations crj on 
the Jeans scale (Meiksin, Bryan & Machacek 2001). Models 
1-3 were run previously to the WMAP results and were con- 
structed to be consistent with existing large-scale structure, 
cluster abundance, and CMB constraints at the time. Model 
4 was later added explicitly to be consistent with WMAP 
results as well (Verde et al. 2003, Hinshaw et al. 2003, Kogut 
et al. 2003). Although Models 1-3 are not the best matches 
to the recent WMAP data, they arc close to models that 
are good matches on the scales relevant to the Lya forest. 
We provide in Tabic 2 a list of similar models with their 
agreement to the WMAP TT data. Model C is the concor- 
dance model of Ostriker & Steinhart (1995). A description 
of the parallel PM code used is given in Paper I. In each 
case we used 512^ particles and a 1024^ force mesh, in a cu- 
bic box with (comoving) side length 25/i~^Mpc, except for 
Model C, for which the comoving box size was 30/i-^Mpc. 
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Qm ^bh^ h n T as Tj 2 {exp(-T)) r_i2 ro A'o JVb, 



0.30 0.70 0.020 0.70 0.95 0.10 0.912 1.31 1020.5 

0.30 0.70 0.020 0.70 0.95 0.15 0.958 1.38 1030.7 

0.30 0.70 0.020 0.70 1.00 0.10 0.954 1.46 1040.8 

0.30 0.70 0.020 0.70 1.00 0.15 1.000 1.54 1041.0 

0.30 0.70 0.022 0.70 0.95 0.15 0.920 1.32 990.1 

0.30 0.70 0.022 0.70 1.00 0.10 0.917 1.40 992.7 

0.30 0.70 0.022 0.70 1.00 0.15 0.961 1.46 990.7 

0.30 0.70 0.024 0.70 1.00 0.10 0.881 1.33 984.0 

0.30 0.70 0.024 0.70 1.00 0.15 0.921 1.39 979.7 

0.30 0.70 0.020 0.75 0.95 0.10 1.068 1.59 1212.1 

0.35 0.65 0.020 0.65 0.95 0.10 0.901 1.32 1001.2 

0.35 0.65 0.020 0.65 0.95 0.15 0.942 1.38 1008.2 

0.35 0.65 0.020 0.65 1.00 0.10 0.940 1.47 1042.3 

0.35 0.05 0.021) 0.()5 1.00 0.15 0.<),S2 1.51 lO.SO.l 

0.35 0.65 0.022 0.65 0.95 0.15 0.907 1.32 986.2 

0.35 0.65 0.022 0.65 1.00 0.10 0.903 1.40 1013.8 

0.35 0.65 0.022 0.65 1.00 0.15 0.946 1.47 1008.1 



0.35 0.65 0.024 0.70 0.95 0.10 0.987 1.49 1022.0 
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Table 3. Some of the output times for the simulations. For 
each redshift the mean Lya flux, (exp(— r)), the required 
H I photoionization rate, in units of 10~^^s~'^, to recover the 
mean flux assuming a homogeneous UV background, the attenu- 
ation length ro, in comoving fe~^Mpc, and the mean number of 
sources per attenuation volume, iVo, are listed for Model 4. The 
last column, A^boxi shows the number of sources in the simulation 
volume of side 25 /i~^Mpc (comoving). The two values at 2 = 4.0 
correspond to two disparate estimates of the mean Lya flux. 

Doppler parameter, where mn is the mass of a hydrogen 
atom. The flux at velocity u is exp(— r). The integration 
variable x indicates the distance along the box in terms of 
the expansion velocity across the box. In evaluating Eq. (1) 
we assume a power-law relation between the gas density and 
temperature 

(2) 

In practice, To and 7 will depend on the reionization history 
of the IGM, but we shall fix To = 2 x lO^K and 7 = 1.5 
throughout as described in Paper I. 

In terms of the baryon density parameter fi^, hydrogen 
baryonic mass fraction X, the metagalactic photoionization 
rate r-12, in units of 10^^^ s~^, and the number /e of elec- 
trons per hydrogen nucleus, the function A{x) is given by: 

A{x) = 18.0/e (^) V-_\^{x)L^^^{l + zf, (3) 

where r_i2(a;) is evaluated at position x, Lbox is in Mpc and 
6 in eq. (1) is in km s~^ . We normalise the spectra according 
to measurements of the mean absorbed flux (exp(— r)) as 
reported in Fan et al. (2002) (their Figure 1) for 2: > 4. We 
discuss the values adopted at other redshifts in Appendix 
B. Values for {exp(— r)) at each of the output times of the 
simulation are provided in Table 3. 

In order to compute the flux power spectra, it is nec- 
essary to filter the simulation data to reduce the efi'ects of 
particle noise at the grid level. The minimal amount of fil- 
tering we consider is Gaussian smoothing on the scale of 
a mesh cell. While this is appropriate for computing the 
power spectrum of the dark matter, it neglects any 'natural' 
filtering arising from pressure forces acting on the baryons. 
The most accurate means of predicting the power spectrum 
is by performing full hydrodynamical simulations at both 
high spatial resolution and in a box of sufficient size to en- 
sure convergence. Such calculations, however, are still pro- 
hibitively expensive to repeat for a range of models, and 
also because the origin of the line broadening is not fully 
understood (Theuns et al. 1999; Bryan & Machacek 2000; 



Table 2. Parameters for cosmological models consistent with 
WMAP data. Q,m is the total mass density parameter, n„ the 
vacuum energy density parameter, fJ;, the baryonic mass fraction, 
h = ifo/100kms~^ Mpc~^ , where Ho is the Hubble constant at 
2 = 0, n the slope of the primordial density perturbation power 
spectrum, r is the Thomson optical depth to the last scattering 
surface, erg the fluctuation normalization at z = in a sphere of 
radius Hh^^ Mpc, ctj the fluctuation normalization at 2 = 3 fil- 
tered on the Jeans scale for gas at T = 2 x 10** K, and the values 
of are from flts to the WMAP TT data. 



Convergence tests, described in Appendix A, suggest that 
these scales are adequate for our purposes. 

Given a set of final particle positions and velocities, 
we compute the spectra as follows. First the density and 
density-weighted line-of-sight velocity are computed on a 
grid (using GIG interpolation) and smoothed using FFT 
techniques. We try three different smoothing schemes, as 
described below and in Appendix A. The default scheme 
simply smooths with a Gaussian of width one grid cell. This 
forms the fundamental data set. A regular grid of 32 x 32 
sightlines is drawn through the box, parallel to the box sides. 
Along each sightline we integrate (in real space) to find t{u) 
at a given velocity u. Specifically we define 



t(u) = / dx A{x 



p{x) 



-\ 2 



T(x)-°-Vie-("-"°^ (1) 



where uo = xaHLhax + wios and b = •v/2fcsT/mH is the 



0.35 0.65 0.024 0.65 1.00 0.10 0.866 1.34 1025.5 
0.35 0.65 0.024 0.65 1.00 0.15 0.907 1.40 1017.6 



To ^ 
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Meiksin et al. 2001), requiring additional modelling which 
affects both the shape and amplitude of the predicted flux 
power spectrum. As we are primarily interested in the rela- 
tive effect of fluctuations in the UV background, we do not 
consider an exhaustive range of smoothing procedures, but 
only consider a few straightforward ones to judge how im- 
portant smoothing is to our final results. Following Gnedin 
& Hui (1998), in addition to the minimal smoothing mod- 
els, we create models with additional smoothing on scales 
of the order of the Jeans length. Specifically we compute a 
'Jeans scale', kj = {3/2y^^c~^H{z)/{l + z), where Cs is the 
sound speed, using the assumed temperature and equation 
of state of the gas and used this in place of our one mesh cell 
filtering above. We consider two different types of smooth- 
ing, a Gaussian of width a = and a form motivated 
by linear perturbation theory, viz. (1 -|- [k/kj]^)~^ . We em- 
phasize again that all of these approaches, like the pseudo- 
hydrodynamical approach employed by McDonald (2003), 
are ad hoc and can only be justified by comparison to true 
hydrodynamic simulations, but modify the power mostly on 
scales where thermal broadening is also very important. 



3 UV BACKGROUND 

3.1 UV intergalactic attenuation 

For the uniform background case, normalising the simula- 
tion results by the mean Lya flux also fixes the attenuation 
properties of the IGM since the ionization structure of the 
IGM has been fixed by the normalisation. This was exploited 
in Paper I to place constraints on the attenuation length ro 
and mean emissivity required to recover the measured values 
of the mean Lya flux for « > 4. At lower redshifts, the sim- 
ulation results arc no longer sufficient for determining the 
Lyman limit optical depth of the IGM. We find that the sim- 
ulations fail to reproduce the measured number of Lyman 
Limit Systems (LLSs) by a factor of several. A similar con- 
clusion was reached by Gardner et al. (1997) on the basis of 
lower resolution simulations than those we present. There, 
the authors attempted to correct for the missing LLSs by es- 
timating the contribution from unresolved halos on the basis 
of Press-Schechter theory. They found that they still fell well 
short of the observed number. Our simulations resolve ha- 
los down to masses of about 10** — 10® M©, corresponding 
to circular velocities smaller than 20 — 25kms~^ at z = 3. 
The baryons in smaller mass halos would be heated to tem- 
peratures too high to remain bound after photoionization 
(Meiksin 1994), so that smaller mass halos would seem un- 
likely candidates for LLSs. We again find too few LLSs. This 
suggests that the origin of these systems may involve more 
complex hydrodynamical processes than simple collapse into 
dark-matter halos, such as Jeans fragmentation in the halos 
(Meiksin 1994), tidal streaming or galactic winds. 

In order to estimate the attenuation length, we adopt 
the measured abundance of LLSs from Stengler-Larrea 
et al. (1995) for tl > 1, dN/dz = No{l + z)~' , with 
A^o = 0.25;° i? and 7 = 1.50 ± 0.39 for 0.3 < z < 4. The 
contribution of the LLSs to the effective optical depth of 
the IGM at the Lyman edge is given by (Zuo 1992b) 




Here, the optical depth corresponds to the path between the 
source at Zcm and the position of observation at Zobs- Be- 
cause of the difficulty in measuring the high optical depths 
involved, the actual distribution of individual Lyman limit 
optical depths is less certain. Adopting a power law distribu- 
tion d^/dzdTL = A(l + z)^T^'^ , and assuming a flat universe 
so that dz/dvp w {HoQ^J^ / c){l + z)^^'^ where rp is proper 
length, we obtain 

^ « ^M^No{l + zf'+- [1 - + r(2 - A 1)] , (5) 

where T(a, x) is the incomplete F— function. Adopting Om = 
0.3, h = 0.7, A^o = 0.25, 7 = 1.50, we obtain dr^s/drp « 
3.2 X 10"'^T(1 + zf, where T = 0.85 for /? = 2 and 0.91 for 
P = 1.5. In the limit of an infinitely steep distribution, so 
that all systems have tl = 1, T = 0.63. 

To estimate the total Lyman limit optical depth 
through the IGM, we add the LLS contribution to the values 
we find from a given simulation for 2: < 4. Current measure- 
ments do not permit an estimate of the LLS contribution at 
higher redshifts. Given the uncertain origin of these systems, 
we conservatively do not include them in our estimates for 
z > A. Doing so would result in a decrease of the attenu- 
ation length. We also assume the proper source emissivity 
does not evolve in our estimate of the attenuation length. 
Over the range of measured optical depths, this assumption 
only weakly affects our results. 

In Fig.l, we show the redshift evolution of dr^g/drp 
from the simulation results for Model 4, along with the 
range in the LLS contribution based on observations. We 
also show the evolution of the attenuation length, defined 
by ro — {dT^g/drp)~^ , and compare our estimate with pre- 
vious estimates based on direct Lya forest line counts from 
Meiksin & Madau (1993) and Haardt & Madau (1996). 
We have included the LLS contribution to the attenuation 
lengths for z < 4. We also show the effect including the 
LLS contribution, if extrapolated to 2; > 5, would have 
on ro by the open points at z = 5 and 5.5. Our esti- 
mates for the attenuation length lie between the medium 
and low attenuation models of Meiksin & Madau (1993), and 
lie systematically above the estimate of Haardt & Madau 
(1996). The redshift dependence of the (proper) attenua- 
tion lengths we find is fit by ro ^ 6800 h'^ Mpc (1 + z)^^'® 
over 2.75 < z < 5.5 to ~ 20% accuracy. The attenuation 
lengths allowing for an extrapolation of the LLS contri- 
bution to z > 4 are fit over the same redshift range by 
ro ~ 1.7 X lO"* h-'^ Mpc (1 + z)-^-'^ to ~ 10% accuracy. We 
caution against extending the fits outside the specified red- 
shift range. In particular, we find the attenuation length de- 
creases precipitously between z = 5.5 and 6, a consequence 
of the observed sharp rise in the Lya optical depth at these 
redshifts (Fan et al. 2002). 



3.2 Source contribution 

The QSO luminosity function has a slope at the bright end 
varying as (t>{L) oc i^'^i , where 0i = 3.41 according to the 
2dF QSO sample for z < 2.3 (Boyle et al. 2000), and /3i = 
j2.58 according to the Sloan QSO sample for z > 3.5 (Fan 
et al. 2001), although /3i = 3.2 is still permitted by the data. 
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Figure 1 . The (proper) gradient drj^ / drp of the Lyman limit 
optical depth (upper panel) and the (proper) attenuation length 

(lower panel) for Model 4. The points are estimated from the 
simulation results using the observational constraints on the mean 
Lya flux. The contribution from Lyman Limit Systems is included 
in the values of tiie attenuation length for z < 4. The effect Ly- 
man Limit Systems would have on the attenuation length at z = 5 
and 5.5 is shown by the open points. The open square at 2 = 4 
corresponds to an alternative lower estimate for the average ab- 
sorbed flux (see text). The two solid curves in the upper panel 
show the upper and lower limits for the measured contribution 
from Lyman Limit Systems, with the dashed extensions extrapo- 
lations based on the lower redshift measurements. The two dashed 
lines in the lower panel arc based on the low and high attenuation 
models of Meiksin & Madau (1993), and the solid line is based 
on their nicdiuni attenuation model. The dotted line is based on 
the attenuation estimate of Haardt & Madau (1996). 



We use the procedure in Paper I for adjusting the parameters 
of the Boyle, Shanks & Peterson (1988) luminosity function. 



[{L/L-r^ + iL/L-r^y 



(6) 



viz., by fixing /32 = 1.58, in accordance with Boyle 
et al. (2000), and matching the SDSS counts for Muso < 
-25.7 for z < 6 and for M1450 < -27.1 at z = 6 (and al- 
lowing for differences in the assumed values of the Hubble 
constant). We show the predicted contribution of QSOs to 
r-12 in Fig. 2 for /3i = 3.2, and compare these with the re- 
quired rates from Table 3. We have followed the procedure 
of Paper I for estimating r_i2 from the QSO cmissivity, 
including hardening of the metagalactic radiation field on 
fihcring through the IGM (Haardt & Madau 1996). Specif- 
ically, for z < 3 we assume a metagalactic spectral index 
ctMG = 0, while for 3.89 < « < 5.5 we take omg = —1.3, ex- 
cept at z = 5.0 (qmg ~ —0.5) and z = 6 (qmg ~ 0). These 
values were chosen to recover the required photoionization 
rates within 2cr, though usually the addition of some diffuse 
radiation is required. The value of omg is uncertain and can 



Figure 2. A comparison between the estimated H i metagalactic 
photoionization rate r_i2 from QSO sources and the rate re- 
quired to match the measured mean Lyo flux, as estimated for 
Model 4. The filled data points with error bars are from Table 
3. The open data point at z = 4 is for an alternative lower es- 
timate of the mean Lya flux (see text). The estimated rates for 
QSO sources, shown by the filled triangles, assume a QSO lumi- 
nosity function with a steep end slope of /3i = 3.2. The inverted 
filled triangles include a 40% boost, nearly the maximum possi- 
ble contribution from diffuse recombination radiation emitted by 
the IGM. The open triangles and inverted triangles are analogous 
estimates assuming /3i = 3.41. The filled square at 2 = 4 is for 
/3i = 3.2 and the alternative lower mean Lya flux considered. 
The open square includes a 40% boost in this estimate due to dif- 
fuse radiation from the IGM. Also shown are the estimates from 
Haardt & Madau (1996) based on earlier QSO source counts and 
assuming an intrinsic QSO spectral index shortward of the Ly- 
man edge of ag = 1.5 (solid line) and ag = 1.8 (dashed line), 
and which include the diffuse recombination radiation from the 
IGM. 



only be computed fully self-consistently by solving the de- 
tailed radiative transfer problem through the high redshift 
IGM, but we do not expect it to deviate much from the 
values we use. In any case, the spectral index only weakly 
affects the overall ionization rate after integrating over the 
photoelectric cross-section. Significant diffuse radiation re- 
sulting from re-radiation by the IGM through radiative re- 
combinations is expected at these high redshifts (Meiksin & 
Madau 1993; Haardt & Madau 1996), with as much as a 
{a A — aB)/ctA ^ 4:0 — 50% boost in the total metagalactic 
photoionization rate over the direct contribution from the 
QSO sources. Wo show the enhanced rates for a 40% boost 
in Fig. 2. 

The QSO contribution to the metagalactic photoion- 
ization rate lies systematically below the required values at 
z = 3 — 4. Given the uncertainties in the required rates, as 
well as in the QSO counts at these redshifts, it is unclear 
how seriously the discrepancy should be taken. Based on 
a different estimate of the QSO counts, Haardt & Madau 
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Figure 3. The auto-correlation function in the UV background 
predicted from QSO sources with /3i = 2.58 (upper panel) and 
/3 = 3.20 (lower panel). The solid curves, from top to bottom, 
correspond to the redshifts z = 6, 5.5, 5, 4, 3.89, 3, and 2.75. 
The z = 4 solid curve is computed for a (comoving) attenuation 
length of ro = 104/i~^ Mpc and nearly coincides with the z = 3.89 
result. Also shown is the result at z = 4 assuming ro = 80h~^ Mpc 
(dashed line), which nearly coincides with the 2 = 5 result. 



(1996) compute a photoionization rate that well matches 
the required rate at z « 3 (see Fig. 2), although their value 
hes much below the required rate hy z = 4. (A required 
slower rate of decline a.t z ^ 3 than predicted by Haardt & 
Madau was also found by Meiksin, Bryan & Machacek 2001.) 
A boost in the IGM temperature at these redshifts, as ex- 
pected for late He ii reionization (Madau & Meiksin 1994; 
Reimers et al. 1997; Meiksin, Bryan & Machacek 2001), 
would also help ease the discrepancy (because of the result- 
ing reduced hydrogen recombination rate). Finally, if the 
luminosity from QSO sources is beamed, the estimated UV 
contribution would increase like the ratio of the UV beam 
solid angle to the optical beam solid angle, since the QSO 
counts are based on restframe optical or near UV (longward 
of the Lyman edge) surveys. (In principle, this could also 
result in a decrease in the estimated QSO contribution to 
the UV background if the ratio were less than unity.) None 
the less, it is interesting in this context that Steidel, Pet- 
tini & Adelberger (2001) report a large Lyman-continuum 
emissivity from Lyman-break galaxies at z ~ 3.4. Adjust- 

1/2 

ing their value down by the factor Qj^j for Model 4, the 
predicted photoionization rate at z = 3, assuming a non- 
evolving proper emissivity, is r_i2 = 4.4, well in excess of 
the required value in Table 3. Within the context of our 
models, such a high rate would so reduce the neutral frac- 
tion of the H I that the IGM would become significantly 
more transparent to Lya photons than measured. To match 
our estimated 3cr upper limit on the photoionization rate 
would require the Lyman-continuum emitting phase of the 



Figure 4. A simulated spectrum from a 'typical' sightline at z=5 
from Model 1. The upper panel shows the spectrum with a uni- 
form ionizing background (A(x) =const, dashed) and where the 
UV background comes from a single source situated at the center 
of the box (dotted). The lower panel shows 0.1 x p/p (dotted) and 
A{x)/A (dashed) for the single source case. 

Lyman-break galaxy population to last only over a redshift 
interval Az ~ 0.1. Thus the required excess photoioniza- 
tion above the QSO contribution could be accounted for by 
the Lyman-break galaxies, but only if they are a bursting 
population at these redshifts. Based on population synthe- 
sis models of the galaxies, Ferguson, Dickinson & Papovich 
(2002) argue this indeed may be the case. 

3.3 UV background fluctuations 

We compute the auto-correlation function of the QSO con- 
tribution to the UV background, using the formalism of Zuo 
(1992b) (also see Paper I), both for /?i = 2.58 and 3.2. The 
results are shown in Fig. 3. Formally the correlation func- 
tion diverges at zero separation, since the variance in the 
radiation field diverges (Zuo 1992a; Paper I). The typical 
width of an absorption feature is 25kms~^ , so when the 
correlations become small on this scale, we expect them to 
have a negligible influence on the flux auto-correlations or 
power spectra. We find that for 2 > 3, the correlations still 
exceed unity on this scale, so that we may expect the effect 
of UV background fiuctuations to imprint a signature on the 
fiux statistics of the Lya forest. We note that for /3i = 2.58, 
additional sources (like galaxies) must be introduced to re- 
cover the required H i photoionization rate (Paper I), but 
that these sources would only dilute the correlations by a 
factor of 2-3. 

The effect of a discrete source on the Lya flux is shown 
in Fig. 4. Here we artificially suppose that the entire UV 
background flux comes from a single source located in the 
center of the simulation box. For a line-of-sight which passes 
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slightly oflF center the A{x) and spectrum are as shown, along 
with the spectrum that would have resulted along this sight- 
line if the UV background were completely uniform. 

To model the fluctuations in the UV background from a 
'realistic' population of discrete sources we resort to a Monte 
Carlo procedure. We create random source lists, with lumi- 
nosities drawn from the QSO luminosity function above. A 
random number of QSOs is chosen consistent with the lu- 
minosity density of QSOs with absolute magnitude at rest- 
frame 1450A Mi45o < -18.5. The luminosity of each QSO 
is drawn at random from the luminosity function, and the 
sources are distributed at random positions in a volume 
larger than the size of the simulation box. We choose the 
box size to be sufficiently large to include the full effects of 
attenuation but not so large that redshift or evolution ef- 
fects become important, so that the box is not laxger than 
the efi^ective size over which sources contribute to the local 
UV intensity (no larger than lu in Paper I). 

The frequency specific flux at the Lyman edge is then 
calculated at each position on the simulation grid from 



i/ro 



(7) 



where the sum is over the sources, the specific luminosity 

of source i at the Lyman edge is Li, and is the distance 
of the ith source from the fiducial point in the siirmlation 
volume. To convert from the QSO luminosity at 1450A to 
the Lyman edge, we adopt a QSO spectral shape of the 



form ~ from 1050 - 1450A and ~ i/" 



from 



912 - 1050A (Zheng ct al. 1997). Since our simulation box, 
and the relevant attcrniation lengths, is small we do not 
keep track of source lifetimes or lightcoue effects (see Croft 
2003). To compute the total rate r_i2 of ionizing photons 
(in units of 10~^^s~'^), from fhix), we adjust the spectral 
index of the metagalactic background to reproduce the tar- 
get values of (exp(— r)) in Table 3. As in the uniform UV 
background case, for z > 3.89 the required values are gen- 
erally negative (—1.3 ^ q;mg ~ — 0.5), consistent with the 
expected hardening of the radiation field near the Lyman 
edge on being filtered through the IGM. Other sources, like 
gala^xies, may of course contribute as well. Adjusting omg 
to reproduce the required ionization rate provides an upper 
limit to the possible effects of UV background correlations 
resulting from attenuation; the introduction of additional 
sources would tend to dilute these effects. We also consider 
cases with a diffuse radiation field, £is is expected to arise 
from radiative recombinations within the IGM. 

Once we have r_i2 throughout the simulation volume, 
we compute A{x). Wo ignore the fluctuations in the IGM 
temperature that would result, noting that this is another 
possible contribution to the flux power spectrum, but one 
again that requires detailed modelling of the reionization 
process. Since our simulation volume is small, we explicitly 
checked that imposing periodic boundary conditions on the 
flux distribution modified our results by much less than 1 
per cent. 

We run sufficient realizations of this procedure to ensure 
convergence (one to several hundred, depending on redshift), 
keeping track of the UV background and flux correlation 
functions and the flux power spectrum. This enables us to 
compute both the mean spectra and their correlations. 
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Figure 5. (Top panel) The cumulative distributions of the pixel 
flux at 2 = 5 for the four models assuming no UV background 

fluctuations. The models arc labelled as Model 1 (dotted), Model 
2 (short dashed). Model 3 (dot-short dashed), Model 4 (long 
dashed). The distributions group according to the value of crj, 
flattening as crj increases. (Bottom panel) The diff'crenccs be- 
tween the cumulative flux distributions for Models 1-3 and the 
distribution for Model 4 assuming no UV background fluctuations 
at z = 5. (The line types correspond to those in the top panel.) 



4 RESULTS 

4.1 Pixel flux distribution 

The cumulative pixel flux distributions at « = 5 are shown 
in Fig. 5. The slopes of the distributions are determined es- 
sentially by the value of aj, and show a systematic trend 
of flattening with increasing aj. The same trend was found 
by Meiksin et al. (2001) at z = 3, who exploited it to 
place stringent constraints on aj from high spectral reso- 
lution measurements of the Lya forest. These results cor- 
roborate the trend found there. For a single high resolution 
spectrum of typically A''pix = 10 pixels, it is possible to dis- 
tinguish the cumulative distributions only to differences of 
about l/ATp/^^ = 0.01, or somewhat larger allowing for pixel 
flux correlations (Meiksin, Bryan & Machacek 2001), so that 
while Models 1 and 3 could be distinguished from Model 4, 
Model 2 could not on the basis of a single spectrum. 

The effect of fluctuations in the UV background on the 
cumulative flux distributions for Model 4 at z = 4, 5, 5.5 
and 6 is shown in Fig. 6. As found in Paper I, the UV back- 
ground fluctuations produce only a small distortion on the 
pixel flux distribution. The distortions found here are some- 
what larger than found in Paper I, because here we now also 
account for correlations in the UV background. As found in 
Paper I, the effects of the UV background fluctuations are 
generally smaller than or comparable to the uncertainty be- 
tween models with differing values of aj. Only at the small- 
est flux values (< 0.1) would it be possible to distinguish the 
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Figure 6. The differences between tlie cumulative flux distribu- 
tions for Model 4 with and without UV background fluctuations, 

at z = 4 (short-long dashed), 5 (short dashed), 5.5 (long dashed) 
and 6 (solid). Also shown is the difference at 2: = 6 allowing for a 
diffuse component added to the fluctuating UV background (dot- 
short dashed). 

effects of the fluctuations from an uncertainty in the model 
if the model is known with sufficient accuracy. A comparable 
level of uncertainty is introduced by the uncertain equation 
of state of the gas (Paper I), but this may be removed by 
adjusting the temperature to match the measured distribu- 
tions of Doppler parameters or wavelet coefficients (Meiksin 
2000; Meiksin, Bryan & Machacck 2001). 

As in Paper I, we find that the demands on the moan 
photoionization rate increase in the presence of UV back- 
ground fluctuations. Allowing for the correlations in the UV 
background, however, reduces the effect, and it actually in- 
verts hy z = 6. Specifically, for z < 5.5 we find the mean pho- 
toionization rate must increase by 5-10% when the (corre- 
lated) UV background fiuctuations are included in the mod- 
els. The amount of the increase depends on the fractional 
contribution arising from a diffuse radiation field and on the 
rcdshift. At 2; = 6, the required moan photoionization rate is 
instead reduced by 10-20% when the now strongly correlated 
UV background fluctuations are included. 

4.2 Pixel flux power spectrum 

We compute the flux power spectrum at wavenumber k for 
the spectra according to 

PF{k) = {\5ff)/{ff, (8) 

where 5f denotes the Fourier transform of 5/ = / — (/). 
In order to assess the accuracy of our predicted flux power 
spectra, we compare the results for Models 1 and 4 and 
the concordance model with the measurements of McDon- 
ald et al. (2000) &t z = 3.89, 3.00 and 2.41, using spectra 
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Figure 8. The P^'' degree of freedom between the concor- 
dance model and the measured flux power spectrum Pp{k), as a 
function of (comoving) resolution. The upper panel is for z=2.7, 
and the lower for z=3.9. The open symbols use all points with 
k < 0.015 (kms^-'^)^-'^ , and the filled symbols use all points 
with k < 0.032 (kms"i)"l . The results are shown for the 
30 h~^ Mpc (comoving) box simulations, for Gaussian smoothing 
on the mesh scale (triangle), Gaussian Jeans smoothing (inverted 
triangle), and Gaussian linear smoothing (square). Also shown 
is the Gaussian mesh-size smoothing case for the 60 h"^ Mpc 
(comoving) box simulation for k < 0.015 (kms~^)~^ (-|-) and 
k < 0.032 (kms~^)~^ (x), as well as the results using alternative 
equations of state (open and filled diamonds) for the 30 h"'^ Mpc 
high resolution box (see text). To help guide the eye, filled sym- 
bols are connected by solid lines and open symbols by dashed 
lines. 



normalised to the identical mean Lya flux they report. The 
results are shown in Fig. 7. The same models are compared 
with the data of Croft et al. (2002b) at 2.72 (using the 
simulation results at z = 2.75). Since Croft et al. report 
no mean Lya flux measurement for their spectra, we adopt 
the value in Table 3. Wc note that according to cq.(26) of 
Croft et al. , Model 1 is ruled out at the 3.9a level, while 
Model 4 and the concordance models at only the l.Ocr and 
1.9(7 levels, respectively. At z — 2.7, we find Model 1 fits 
about as well as Model 4 and the concordance model. For 
k < 0.015 (kms^^)"^ , the for Models 1, 4 and the con- 
cordance model, compared with the data of Croft et al. , 
are, respectively, 13.8, 16.5 and 18.3. The respective prob- 
abilities P(> x^) for the 8 fc-valuos are 0.087, 0.036 and 
0.019, so that Model 1 is actually preferred. 

In order to understand the differences between our re- 
sults and those of Croft et al. , we have performed a res- 
olution sequence of comparisons between the concordance 
model simulations and the measurements, for a fixed equa- 
tion of state (ic, leaving To and 7 unchanged). In Fig. 8 we 
show how the reduced x^ varies with resolution, for both 
k < 0.015 (kms-^)-^ and k < 0.032 (kms-^)-^ (The simu- 
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Figure 7. The flux power spectrum Pp{k), as a function of restframe velocity wavenumber, for Models 1 (x) and 4 (square) and the 
concordance model (inverted Y), compared with the measurements in high resolution spectra at the indicated redshifts. 



lation values are linearly interpolated to the fc— values of the 
measured Prik) points.) We do the comparison for mesh- 
scale smoothing, Gaussian Jeans smoothing and linear Jeans 
smoothing. A clear trend of diminishing x with increasing 
resolution is apparent in all cases, showing that the lower 



resolution simulations have not yet matched the precision of 
the measured Ppik) values. Only at the highest resolution 
do the simulations begin to converge, giving values of suf- 
ficient precision to make a meaningful comparison with the 
measured values. The decrease in a-t higher resolution is 
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Figure 9. The flux power spectrum PF{k), as a function of restframe velocity wavenumber, for Model 4, with (square) and without (x) 
UV background fluctuations. Also shown is Ppik) at z = 6 allowing for a slight suppression due to diflFuse radiation (inverted Y). The 
results are shown at z = 4, 5, 5.5 and 6. 
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Figure 10. The ratio of tlie flux power spectrum PF{k) with UV background fluctuations to the value without, as a function of restframe 
velocity wavenumber, for Model 4. The results are shown at z = 2.75, 3, 3.89 and 4. 



noteworthy, as it suggests a given model is a better fit to the 
data than would be believed on the basis of lower resolution 
simulations. The larger values found for the 60 Mpc 
box for k < 0.032 (kms~^)~^ suggest that the simulations 



have still not quite converged on these scales, at least at a 
resolution of Ar — 58.6 h~^ kpc, particularly for z < 3. 

The results in Fig. 8 also show that the goodness-of- 
fit is highly sensitive to the smoothing. Adopting Gaus- 
sian Jeans smoothing instead of smoothing on the mesh 
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Model P{> X^) 



1 


55.1 


0.007 


2 


38.7 


0.192 


3 


42.3 


0.106 


4 


40.7 


0.139 


C 


44.9 


0.065 



Table 4. Statistical comparison betwen the model predictions 

and measurements for the flux power spectrum, combining data 
at z = 2.41, 2.72, 3.00 and 3.89. No adjustment to the equation 
of state or mean Lya flux have been made. 



scale results in = 36.6 for k < 0.015 (krris^^)^^ , and 
P{> X^) = 1-4 X 10~^, in considerably poorer agreement 
than the mesh-scale smoothing. The agreement is also much 
poorer than predicted by Croft et al. , although we have 
not attempted to adjust the mean flux or equation of state 
to improve the match. Similarly poorer results arc found 
for linear Jeans smoothing. We also compute the effect of a 
change in the equation of state for the mesh-scale smoothing 
case, adopting the values for To and 7 predicted by the late 
Ife II reionization model of Schaye et al. (2000). The agree- 
ment considerably worsens, as shown in Fig. 8, resulting in 
values comparable to the Jeans smoothing cases. 

Without direct comparisons to full hydrodynamical 
simulations of comparably high resolution, it is unclear how 
much artificial smoothing should be applied to the PM re- 
sults, or, indeed, how accurately any smoothing scheme may 
mimic the true underlying hydrodynamical processes. A cer- 
tain amount of smoothing will always be present in the spec- 
tra resulting from thermal broadening and peculiar velocities 
in the gas, which the PM simulations are able to account for. 
Indeed, there is a certain amount of degeneracy between the 
Jeans-scale broadening and changes in the equation of state. 
The crucial question is whether the pressure forces missing 
in the PM simulations result in significantly less smoothing 
than should be physically present to obtain accurate pre- 
dictions for Ppik). We are unable to address that question 
here, but we have demonstrated that the successful applica- 
tion of PM (or HPM) simulations for discriminating between 
cosrriological models on the basis of the fiux power spectrum 
relies critically on the adopted degree of artificial smoothing. 

We compare our models with the data of McDon- 
ald et al. and Croft et al. (2002b) in Table 4 combin- 
ing the results at a = 2.41, 2.72, 3.00 and 3.89 for k < 
0.015 (kms~^)~^ . The total number of degrees of freedom 
is 32. We have conservatively adopted mesh-scale smooth- 
ing throughout. Within the context of mesh-scale smooth- 
ing, the results we provide are upper bounds to the for 
each model, since we have made no effort to improve the fits 
by making the several adjustments available. Specifically, we 
have not: 1. adjusted fa to improve the fits, adopting instead 
the values given by McDonald et al. , or from Table 3 for 
z = 2.7, 2. adjusted the equation of state parameters To and 
7, 3. accounted for the spread in the model predictions of 
Ppik), or 4. accounted for mode-mode correlations in the 
power either for the simulations or the data. We confine the 
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Figure 11. (Upper panel) The ratio of the flux power spectrum 
for Model 1 (5-pt star). Model 2 (x) and Model 3 (inverted Y) 

to that of Model 4 at z = 5. (Lower panel) The ratio of the flux 
auto-correlation function for Model 1 (dotted). Model 2 (short 
dashed) and Model 3 (dot-short dashed) to that of Model 4 at 
z = 5. 



comparison to < 0.015 (kms~^)~^ , since for this range 
the model predictions are more stable to changes in the box 
size, and the range avoids an anomalous low excursion in 
the measured Ppik) at a = 3 for fe = 0.0179 (kms"^)"^ with 
very small error bar. 

The effect of UV background fiuctuations on the fiux 
power spectrum is shown in Figs. 9 and 10. Since the UV 
background fluctuations tend to somewhat suppress the 
mean flux (boost the mean optical depth) (Paper I), we have 
re-normalised the uniform background case to the identical 
mean flux found for the fluctuating background case at each 
redshift to make a fair comparison. The effect of the fluctu- 
ations are most pronounced on large scales and small scales. 
The small scale peak is on the scale of the absorption line 
widths and is due in part to a resetting of the mean ioniza- 
tion rate necessary to reproduce the same mean Lya flux. 
At intermediate scales the power tends to be somewhat sup- 
pressed. The trend at low k of decreasing power with in- 
creasing k continues even to z = 2.75, although the effect 
is much reduced and essentially vanishes by 2 = 2.41. The 
rise at the lowest k value appears to be real, not an artefact 
of the flnite box size, based on convergence tests applied to 
the concordance model (see Appendix A). At z = 6, we also 
show the small dilution effect of introducing an additional 
diffuse radiation background of 50% of the average direct 
contribution from the QSO sources. 

The ratios of the flux power spectra for Models 1-3 to 
Model 4 are shown in Fig. 11 at z = 5. The differences be- 
tween the models are smaller than the effect induced by the 
UV background fluctuations on the largest scales (smallest 
wavenumbers) . The effect of the fluctuations should be dis- 
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Figure 12. The flux auto-correlation function ^piv) at 2: = 2.7, 
as a function of restframe velocity separation, for Model 4 (solid) 
and the concordance model (dashed), assuming Gaussian smooth- 
ing on the mesh scale. Also shown is the concordance model using 
linear Jeans smoothing (dotted). The models are compared with 
the measurements from Croft et al. (2002b). 



Figure 13. The flux auto-correlation function ^f(v), as a func- 
tion of rest frame velocity separation, for Model 4, with (solid) 
and without (dashed) UV background fluctuations. Also shown 
is ^p{v) at 2; = 6 including UV background fluctuations, but al- 
lowing for a slight suppression due to additional diffuse radiation 
(dot-dashed line). The results are shown aX z = 4, 5, 5.5 and 6. 



cernable in z > 5 spectra of sufficient signal-to-noise ratio 
within the current precision of the measured cosmological 
parameters. 



4.3 Pixel flux auto-correlation function 

We compute the flux auto-correlation function at a velocity 
separation v for the spectra according to 

^Fiv) = {f{v')f{v'-v)}/{ff-l, (9) 

where f{v) denotes the pixel flux at velocity v, and (/) is 
the measured mean flux at the relevant redshift. The results 
are shown for Model 4 and the concordance model in Fig. 12 
at z — 2.75. The data points arc from Croft ct al. (2002b) 
at (z) = 2.72. While the agreement is close at small sepa- 
rations, the model estimates disagree significantly with the 
measured values. The agreement worsens at large separa- 
tions, likely due to lack of numerical convergence in the sim- 
ulation results. We have not attempted to adjust the mean 
Lya flux or the equation of state to improve the match. 

We compare the flux auto-correlation function with and 
without UV background fluctuations in Fig. 13. The largest 
effect is at 2: = 6 and decreases with decreasing redshift. 
The effect almost disappears by z = 4 for velocity separa^ 
tions smaller than 200 km s~^ . We note that the convergence 
of the auto-correlation function is not guaranteed for veloc- 
ity separations exceeding ~ 100 km s^"*^ . At 2: = 6, we also 
show the small dilution effect of introducing an additional 
diffuse radiation background of 50% of the average direct 
contribution from the QSO sources. 



The ratios of the flux auto-correlation function for Mod- 
els 1-3 to Model 4 arc shown in Fig. 11 at z = 5. The differ- 
ences between the models are smaller than the effect induced 
by the UV background fluctuations for small to moderate 
scales (< 100 kms"^ ). 



5 CONCLUSIONS 

We have examined the effect of UV background fluctua- 
tions expected from QSO sources on the power spectrum 
and auto-correlation function of the Lya forest. To do so, 
we first performed several sequences of tests to assess the 
simulation parameters required to reach convergence. We 
investigate the convergence for different amounts of smooth- 
ing meant to mimic hydrodynamical effects, almost all for a 
fixed equation of state. We note that varying the equation 
of state will affect the predicted flux power spectrum, and 
should be included when attempting to constrain cosmolog- 
ical models by comparison with the data. We have not per- 
formed an exhaustive search in smoothing algorithms and 
equations of state to achieve a best match to the data. Our 
emphasis in this paper is instead to explore the impact these 
have on the relative effect UV background fluctuations are 
expected to have on the predicted flux power spectrum and 
auto-correlation function. Our conclusions concerning these 
tests are (where all length scales are comoving and velocity 
scales proper): 

1. To reach 10% convergence in the required 
H I photoionization rate r_i2, a (comoving) spatial resolu- 
tion of 60 kpc is required. 
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2. We were not able to reach convergence in the ID dark 
matter power spectrum to better than 50% on any scale. 
Box sizes exceeding 60 Mpc and/or a spatial resolution 
of better than 30 kpc appear to be required at 2 = 5. 
The convergence requirements arc oven more severe at lower 
redshifts. Fortunately good convergence in the ID matter 
power spectrum is not required for good convergence in the 
predicted flux power spectrum because thermal and Doppler 
broadening of the gas filters out the contribution of the high 
spatial frequency modes to the H i optical depth. 

3. We show that the predicted flux power spectrum 
Ppik) is senstitive to the degree of smoothing. A crucial 
question is how much (or even whether) any additional 
smoothing to the gas variables should be applied beyond 
the inherent smoothing of the spectra due to Doppler and 
peculiar velocity broadening. For the minimal smoothing on 
the scale of a mesh cell, the minimum resolution and box 
size required to achieve 10% convergence on (proper) scales 
k < 0.015 (kms"^)"Ms 30 /i"^ kpc in a 25ft"^Mpc box. 
We were not able to achieve convergence to 10% at larger 
wavcnumbers. Allowing for additional smoothing improves 
the convergence. With the application of Gaussian Jeans 
smoothing, it is possible to achieve convergence on all scales 
A: < 0.1 (kms~^)~^to 10% at 2 = 5 for a spatial resolution 
of at least 60 kpc in a 30 Mpc box. With the appli- 
cation of linear Joans smoothing, we arc able to achieve this 
level of accuracy only for k ^ 0.01 (kms"'^)^^ . The conver- 
gence requirements become more severe at lower redshifts. 

4. We find that introducing Jeans Gaussian or Jeans lin- 
ear smoothing significantly flattens the predicted flux power 
spectrum near its peak. Without full hydrodynamical sim- 
ulations at comparable resolution to our highest resolution 
simulations in a box size of at least 25 Mpc , it is not pos- 
sible to assess whether or not the additional Jeans smooth- 
ing, while it improves the convergence properties, also im- 
proves the accuracy of the flux power spectrum as predicted 
from pure gravity simulations. 

5. We were unable to obtain convergence in the flux 
auto-correlation function ^F{k) to better than 10% on spa- 
tial separations exceeding ~ 3% of the box size. Compar- 
isons with the data for this statistic appear problematic. 

6. A criterion of the level of simulation convergence 
needed to make a meaningful comparison with observations 
is the convergence in between the simulation predictions 
and the measured values. We find that for the com- 
parison between the predicted and measured Ppik) contin- 
ues to decrease with increasing resolution. For our minimal 
mesh-scale smoothing case, we find differences in the re- 
duced 2 f5i 3 — 4 of ~ 50% on going from a reso- 
lution of ~ 60/i~^kpc to ~ 30/i~^kpc in a 30/i~^Mpc 
box, and differences of ~ 20% on increasing the box size to 
60/i~^Mpc at a fixed spatial resolution of ~ 60/i~^kpc. 
Allowing for Jeans smoothing does not greatly improve the 
convergence in the reduced confined to wavenumbers 
k < 0.015 (km s^^)^^ . In particular, we find that while the 
concordance model would be very strongly rejected on the 
basis of moderate resolution simulations using the minimal 
amount of mesh-scale smoothing, the model is rejected at 
only the ^ 2a level in our highest resolution simulations. 
By contrast, when Jeans smoothing is applied, the model 
is very strongly rejected at all resolutions. We have not ex- 



plored variations in the equation of state that could poten- 
tially improve (or degrade) the agreement. 

7. We find significant degeneracy between the smooth- 
ing and the equation of state of the baryons. It is unclear 
how the degeneracy may be removed without performing full 
hydrodynamic simulations to infer the equation of state of 
the IGM by comparing with the measured small scale fluc- 
tuations of the Lya forest on the scale of the absorption fea- 
tures. A possibility is to use the measured Doppler parame- 
ters to constrain the equation of state (Schayc ct al. 2000), 
but without high resolution simulations in box sizes of at 
least ~ 25 Mpc , it is unclear whether or not the param- 
eters deduced in this way are the optimal ones to use for pre- 
dicting the flux power spectrum. One complication is the ev- 
ident need for additional broadening at z < 3.5 beyond sim- 
ple photoionization heating (Theuns et al. 1999; Bryan & 
Machacek 2000; Meiksin, Bryan & Machacek 2001). The dis- 
crepancy in the predicted and measured line widths suggests 
the underlying physics is still not fully understood, partic- 
ularly for the moderate to low optical depth lines (Meiksin, 
Bryan & Machacek 2001). 

The simulations are normalised on the basis of the av- 
erage Lya flux transmitted through the IGM. We find that 
there is significantly more accord in the measurements of 
this parameter than previously recognized once existing con- 
straints are interpreted on a consistent statistical basis. 

We assess the possibility of using measurements of 
the fiux power spectrum for inferring the contribution of 
QSO sources to the metagalactic UV ionizing background. 
Because of the sparsity of QSOs relative to alternative 
sources, like galaxies, the UV background contribution from 
QSOs will exhibit large, spatially correlated fluctuations. 
We model the fluctuations using Monte Carlo realisations 
of randomly distributed QSO sources based on QSO counts 
from the 2dF and SDSS QSO surveys. We estimate the at- 
tenuation lengths through the IGM self-consistently from 
our simulations, requiring a background ionization rate that 
recovers the measured mean Lya flux. Because the simular- 
tions are not able to recover the measured number of Lyman 
Limit Systems for 2 < 4, we add their contribution to the 
attenuation by hand. 

The resulting fluctuations increase the required pho- 
toionization rates by 5-10% over the uniform background 
case for z < 5.5 in order to match the measured mean 
Lya fluxes. This is a somewhat smaller effect than was found 
in Paper I, where only the UV background fluctuations were 
accounted for and not their spatial correlations. By 2 = 6, 
we find the trend reverses in the presence of the now very 
strong UV background correlations: the required photoion- 
ization rate is reduced by 10-20% compared with the uniform 
background case. By contrast, we find that the effect on the 
cumulative flux distribution increases when the correlations 
in the UV background fluctuations are included. Previously 
we found the fluctuations distorted the flux distributions by 
~ 0.2% for 2 > 5. Allowing for the correlations increases the 
effect to 1-1.5%, rendering them somewhat easier to detect. 

We find that the magnitudes of the relative effects on the 
2-pt flux statistics by fluctuations in the UV background are 
much better converged than the absolute predictions. The 
ratios of the predicted flux power spectra with and without 
UV background fluctuations are converged to better than 
1% on scales k < 0.2 (kms~^)~^ in our simulations, and are 
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in agreement for Gaussian rncsh, Jeans Gaussian and Jeans 
linear smoothings at this level over these scales. The ratio 
of the predicted auto-correlation function with and without 
UV background fluctuations, however, is converged to better 
than 10% only for velocity separations ^ 200 km s^^ . 

We also find that the UV background fluctuations 
boost the flux power spectrum at laxge scales {k < 
0.01 (kms^^)"^ ) and small scales {k ~ 0.2-0.4 (kms"^)"^ , 
corresponding to the widths of individual absorption fea- 
tures), and suppress the power at intermediate scales. For 
z < 4, the efi^ects are small, at the few percent level. At 
higher redshifts, however, the large scale effect grows due to 
the increasing strength of the UV background correlations, 
resulting in a boost in power of over 50% by 2 > 5. By 
2 = 6, the background fluctuations produce an upturn in 
kPF{k)/-iv towards low k, providing distinctive evidence for 
the presence of UV background fluctuations, as such an up- 
turn would be unphysical from large-scale structure effects 
alone. Measuring such an effect at the current level of flux 
detection is problematic because of the very high Lya optical 
depths at these redshifts, but it may become feasible with 
future Extremely Large Telescope technology, or possibly by 
searching for the signal at higher orders in the Lyman series 
where the optical depths are reduced. 




k (h/Mpc) 

Figure Al. The (real space) ID matter power spectrum in the 
concordance cosmology at 2 = 5 from our resolution sequence 
with 128^, 256'^ and 512=^ particles and force meshes of 256^, 512^ 
and 1024^. The solid line shows the predicted ID spectrum, using 
the Peacock & Dodds (1996) fitting formulae for the 3D matter 
power spectrum, Eq. Al, and a Gaussian filtering of 3 mesh cells 
(FWHM). 
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APPENDIX A: CONVERGENCE TESTS 

Even including only gravitational physics, simulations of the 
Lya forest are restricted by two size constraints. High spa- 
tial resolution is required to resolve the structure of the ab- 
sorbers, while a large box size is required both to capture 
the large-scale power which will affect the velocity widths of 
the absorption lines as well as to obtain a fair sample of the 
universe. Although a small box may be adequate to resolve 
the line structure, if the box is too small the fluctuations 
corresponding to the scale of the box will become nonlin- 
ear, in which case the results are no longer representative of 
the cosmological model simulated. These numerical effects 
have been investigated in different contexts by a number of 
authors (sec references in McDonald 2003 for example). 

To investigate these effects of numerical resolution fur- 
ther we chose a relatively large box, 30ft~^Mpc (comoving) 
on a side, and simulated the concordance cosmology using 
128^, 256^ or 512^ particles (corresponding to 256*, 512* 
and 1024* mesh cells). At the high redshifts of interest to 
us the fundamental mode of the box is safely in the lin- 
ear regime, suggesting that the simulation represents a fair 
sample of the universe. To isolate the effects of numerical 
resolution from sample variance, we sampled the initial con- 
ditions starting at low-fc and moving 'outwards' to higher 
k so that the initial phases and amplitudes of the common 
modes were identical from simulation to simulation. 

The 3D matter power spectrum (not shown) exhibits 
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Figure A2. The ID flux power spectrum in the concordance 
cosmology at z = 5 from our resolution sequence with 128^, 256^ 
and 5123 particles and force meshes of 256^, 512^ and 1024^. The 

box size in each run is 30 h^^Mpc. These results indicate that our 
highest resolution runs are converged at better than the 10% level 
(in power) on scales k < 0.015 (kms~^)~^ . 



the expected behaviour, with the higher resolution simula- 
tions having progressively more small-scale power. At 2 = 5 
the fundamental mode is safely in the linear regime, with 
A^(fcfund) — 0.04, and the mesh scale is non-linear with 
A^(fcniosh) > 10. The fitting formula of Peacock & Dodds 
(1996) underestimates the uon- linear power by 0(50) per 
cent. Fig. Al shows the ID matter power spectrum, in real 
space. This is an integral of the 3D spectrum over wavenum- 
bers greater than the ID wavenumber under consideration: 
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Figure A3. The ID flux power spectrum in the concordance cosmology aX z = 5 from a resolution sequence with 128^, 256^ and 512^ 
particles and force meshes of 256^, 512^ and 1024^. The box size in each run is 30/i~^Mpc. (Left) The density and velocity fields are 
Gaussian smoothed by the Jeans scale. (Right) The density and velocity fields are smoothed by (1 + [fe/fej]^)~^. 
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Figure A4. The ID flux power spectrum in the concordance 
cosmology at z = 3 from a resolution sequence with 128^, 256^ 
and 512^ particles and force meshes of 256^, 512^ and 1024^. The 
box size in each run is 60 h~^Mpc. 
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We find that the numerical results are relatively well fit by 
the Peacock & Dodds (1996) formalism, if we filter the 3D 
power spectrum, A'^jy{k), with a Gaussian filter of FWHM 
3.2 mesh cells (plotted as the solid line for the highest reso- 
lution run) to approximate the finite force resolution of the 
PM code. Note that missing high-fc power in the lower res- 
olution simulations biases the ID power spectrum low over 
the whole fc-range plotted since much of the long-wavelength 
power comes from aliasing of 3D short wavelength power. 
We show in Fig. A2 the power spectra of the flux produced 
by these 3 runs using a Gaussian smoothing of one mesh 
cell to define the density and velocity fields. (The flux nor- 
malisation is slightly different than used in the text.) From 
these results we infer that, for the relatively high redshifts of 
interest, a 30/!-~^Mpc box with 512^ particles and a 1024'' 
force mesh produces spectra converged to better than the 
10 per cent level. 

Similar conclusions were found by McDonald (2003), 
though our convergence is not quite at the level quoted 
therein. Part of the reason for this is the treatment of small- 
scale power in the simulations. To investigate this we modi- 
fied our smoothing procedure following Gnedin & Hui (1998) 
to mimic the effects of thermal pressure on the baryonic 
component, using either Gaussian or linear smoothing as 
described in Section 2 above. 

We show in Fig. A3 that the Gaussian smoothing 
method, used by Zaldarriaga et al. (2001), improves the con- 
vergence. This is to be expected since the Jeans scale here is 
close to the mesh size of our intermediate resolution run. We 
found that pure linear theory smoothing introduced numer- 
ical artifacts unless it was supplemented with our standard 
single mesh cell Gaussian smoothing, so we show in Fig. A3 
the results of this hybrid. For the smaller runs the Gaussian 
smoothing dominates, but at higher resolution the linear 
Jeans smoothing is suppressing small-scale power as well. 



Figure A5. The ratios of the (upper panel) ID flux power spec- 
trum and (lower panel) flux auto-correlation function, with and 
without UV background fluctuations, in the concordance cosmol- 
ogy at 2 = 4 from a resolution sequence with comoving grid 
size Ar = 58.6h~^kpc, in boxes of comoving sides 30/i~^Mpc 
(x) and 60 Mpc (inverted Y), and with a comoving grid size 
Ar = 29.3 /i^^ kpc in a box of comoving side 30h~^ Mpc (5-pt 
star). 



The convergence becomes fractionally better if we in- 
crease To or decrease 7. It becomes only slightly worse if 
we lower To to lO'^K. A further sequence reducing the box 
size (not shown here) suggests that we can safely reduce the 
volume simulated to a box of side 25/i^^Mpc if we are in- 
terested in the Lya forest at high redshift. This makes the 
mesh scale slightly smaller, capturing the features in the for- 
est slightly better, thus we have used 25 /i~^Mpc thoughout. 

Finally, we show in Fig. A4 a convergence study at lower 
redshift {z = 3) in a larger box (60 /i~^Mpc). The power 
spectrum in the larger box follows closely that of the smaller 
box at half the resolution. It thus departs significantly from 
that in the 30 ft~'^Mpc box at fc ~ 0.01 km~^ s for the highest 
resolution runs, slightly earlier than might be expected from 
Fig. A2 due to the increased non-linear scale. 

In Fig. A5, we show a convergence sequence at z = 4 
for models including fluctuations in the UV background. The 
UV background fluctuations are based on the QSO counts 
with l3i = 3.2 described in the text, with an added dif- 
fuse component increasing the ionization rate by 40% to 
reproduce the measured average flux. The rise in Prik) at 
the lowest fc-values appears to be real. The effect of the 
UV background fluctuations on the flux auto-correlation be- 
comes poorly determined for velocity separations exceeding 
300 km s"^ . 

We show in Fig. A6 that the Gaussian and linear 
Jeans smoothing improve the convergence of the flux auto- 
correlation function over the Gaussian single-mesh scale 
smoothing. For 256^ particles on a 512'' force mesh, the 
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Figure A6. The error in the ID flux auto-correlation function in 
the concordance cosmology at z = 4 from (panels a-c) a resolution 

sequence in a box of comoving size 30ft^^Mpc with 128'^ and 
256"' particles and force meshes of 256'^ and 512'^, compared with 
512'^ particles and a force mesh of 1024'^, and from (panel d) a 
box size sequence comparing boxes of comoving size 30/i~^Mpc 
with 60/i~^Mpc for (comoving) force mesh resolutions of Ar = 
117h-iMpc and 58.6 h"! Mpc. 



Jeans smoothings give convergence to better than 5% on 
comoving separations of r < 3/i~^Mpc, or 10% of the box 
size. The singlc-incsh Gaussian smoothing converges to only 
15% on these scales. Also shown is a comparison of the 
single-mess Gaussian smoothing for comoving box sizes of 
30 h~^Mpc and 60 h~^Mpc, at two spatial resolutions, show- 
ing convergence to 10% only for comoving separations of 
r < 1 h~^Mpc, or 3% of the box size. 

While the smoothing schemes introduced above improve 
convergence, this does not guarantee that they coverge to 
the correct power spectra, as would be computed using a 
full hydrodynamical scheme. A direct comparison is still 
prohibitively computationally expensive because of the wide 
range in lengthscales required both to resolve the non-linear 
structures into which the bulk of the IGM collapses and 
to capture sufficient large-scale power to obtain an accu- 
rate computation of the flux power spectrum and flux auto- 
correlation function. The most we can demonstrate is that 
the results are not excessively sensitive to the smoothing 
scheme. 

In Fig. A7, we compare the three smoothing methods 
for the concordance model at z = 4. The power spectra 
for both Jeans smoothings agree to within 1-2%, although 
they are suppressed in power compared with the Gaussian 
single-mesh smoothing by as much as 8% near the peak. 
By contrast, both Gaussian smoothings result in flux auto- 
correlation functions agreeing to 3%, while the linear Jeans 
smoothing results in an excess of correlations of as much 
as 60%. These results leave ambiguous which smoothing 



Figure A7. The ID flux power spectrum (upper panel) and flux 
auto-correlation function (lower panel) in the concordance cos- 
mology at z = 4 for Gaussian smoothing on the mesh size, Jeans 
Gaussian smoothing and Jeans linear smoothing. The box size in 
each run is 30/i~^Mpc. 

method may be most accurate. For simplicity, we retain the 
mesh-scale Gaussian smoothing for the main results in the 
text, but make reference to alternative smoothings. 

The different smoothing methods result in different 
rates of convergence of the metagalactic ionization rate r_i2 
required to match the measured mean Lya flux, as shown 
in Fig. A8. For a force mesh of 1024^ , the different methods 
yield ionization rates agreeing to 2%. This is comparable to 
the difference between meshes of sizes 512"^ and 1024'^ for 
each smoothing, except for the linear Jeans smoothing, for 
which the diflference is 4%. 



APPENDIX B: MEAN FLUX 

The simulations leave unspecified the ionization level of the 
IGM. As a consequence, the flux normalisation of the spec- 
tra is undetermined. We fix the normalisation of the simu- 
lations by adjusting A in eq. 3 to match the measured mean 
Lya flux values as given by Fan et al. (2002) in their Fig- 
ure 1 for 2 > 4, as tabulated in Table 3. For z < 3.89, we 
estimate the mean flux from a set of published Keck HIRES 
spectra, as described below. 

The measurements of the mean flux have proven dif- 
ficult, largely due to uncertainties in the QSO continuum 
across the absorbing region. Deriving error bars on the esti- 
mates is even more difficult because of difficulties in measur- 
ing the error on the continuum estimate, and because of flux 
correlations between pixels. PRS used bootstrap resampling 
to estimate the errors on their fit to the mean Lya optical 
depth fa, obtaining fa = A{1 + zy+^ , with 7 = 2.46 ± 0.37 
and A = 0.0175-0.00567±0.0002. The interpretation of the 
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Figure A8. The convergence of the required metagalactic ion- 
ization rate r_i2, in units of lO"'^^ s"'^, at z = 4 for force mesh 

(comoving) resolutions of Ar = 29.3 kpc, 58.6 kpc and 
117 h^-"^ kpc, for Gaussian smoothing on the mesh size, Jeans 
Gaussian smoothing and Jeans linear smoothing. Results shown 
are for the concordance model. The comoving box size in each 
run is 30 h~^Mpc. 
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Figure Bl. The probability distribution of the mean absorbed 
Lya flux at z = 2.72 according to the estimate of PRS. The long 
vertical tick mark indicates the average mean flux. 
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0.48 ± 0.05 
0.51 ±0.05 
0.71 ± 0.03 

0.76 ±0.03 
0.82 ± 0.02 



0.49 ± 0.01 
0.51 ±0.01 
0.71 ± 0.01 

0.76 ±0.01 
0.82 ±0.01 



0.47 ± 0.03 
0.48 ± 0.02 
0.70 ± 0.02 
0.74 ± 0.04 



Table Bl. Estimates for (exp(— t)) using fits to tq, from Press, 
Rybicki & Schneider (1993) (PRS), Zhang et al. (1997) as fit 
to the Da values of Steidel & Sargent (1987) (SS87), and Kim, 
Cristiani & D'Odorico (2001) (KCD). Also shown are the values 
adopted in this paper. 



error is not completely clear. Seljak, McDonald & Makarov 
(2003) argue that, if the errors from PRS are interpreted 
as Gaussian, then the PRS fit results in both a bias and 
large error estimate for fa. Seljak et al. find at « = 2.72, 
{fa) = 0.29±0.14, where {. . .) denotes an average over Gaus- 
sian distributions for A and 7, while using the mean values 
of A and 7 would instead give fa = 0..35. The error they 
quote, however, is the statistical rms on fa, which carries 
little probabilistic meaning since the resulting probability 
distribution of fa is highly skewed. Moreover, they allowed 
for unphysical values of fa < in their averaging. Truncat- 
ing the distribution to fa > yields an rms of 0.09 instead, 
which we estimate from a Monte Carlo calculation. In terms 
of an equivalent Gaussian la error spread, defined so that 
the probabilities of obtaining a value la above or la below 
the mean are each 0.16, the equivalent la error spread in fa 
is (fa) = 0.30lg;Qg, truncating to fa > 0. 

Relating the mean optical depth to the mean absorbed 
flux and its error is not straightfoward. Instead it is bet- 
ter to start with the distributions for A and 7 and esti- 
mate (exp(— fa)) directly. The resulting probability distri- 
bution for the mean flux at x = 2.72 is shown in Figure Bl. 
The peak of the distribution is at exp(— fa) = 0.711q;q2) 
where the error bars refer to the range in the mean flux en- 
closing 68% of the probability. The moan is (exp(— fa)) — 
0.75 ± 0.07, where the error is the rms, and the distribution 
was truncated to fa > 0. Adopting the la equivalent proba- 
bilistic interpretation of the error, as above, gives instead a 
mean of (exp(— fa)) = 0.75lQ;Qg. Because of the skewed na- 
ture of the probability distribution, the estimate of the error 
for any quantity depending on the mean flux must take into 
account the full probability distribution of the mean flux. 

For 2: < 4, we use the published Keck HIRES spectra 
used by Mciksin, Bryan & Machacek (2001) for testing mod- 
els of the Lya forest. These spectra cover the redshift range 
2.5 < z < 3.7. We estimate the mean flux from each spec- 
trum over a narrow redshift interval {Az = 0.1—0.3) centred 
at each redshift considered, and average the results, adopt- 
ing the error in the mean as our error estimate. We also allow 
for small systematic offsets in the continuum level (typically 
1-4%), as provided by Meiksin et al. for the two best-fltting 
models to the flux distribution (which are ACDM models 
similar to those considered in this paper). We estimate the 
error in the systematic offset from the spread between the 
best-fitting models (1-2% of the continuum level), and add 
this linearly to the error in the mean. (Since the error is sys- 
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tematic rather than statistical, adding the error in quadra- 
ture would be inappropriate.) 

For z — 3.89, wc adopt the value of McDonald 
et al. (2000) of 0.48 ± 0.02. For z = 4.0, we use the best-fit 
of Songaila & Cowie (2002) to the measured redshift trend 
(their oq.[19]), to obtain a value for the mean flux of 0.47. 
We also linearly extrapolate the values of Songaila & Cowie 
(2002) at (z) = 4.34 and 4.07 down to z = 4.0, obtaining 
a mean flux of 0.36 ± 0.03. The difference appears due to a 
single outlying measurement which substantially lowers the 
mean. Wc shall consider both values, noting that the larger 
one better matches the trend in Tabic Bf . 

For comparison, we also show the estimates, using the 
above procedure as applied to the fit to fa of Press et al. , 
based on similar estimates from Kim, Cristiani & D'Odorico 
(200f) with A = 0.0144 - 0.004717 and 7 = 2.43 ± 0.17 
(which we infer from their error estimates for A and 7), and 
the fit by Zhang, Anninos, Norman & Meiksin (1997) to 
the flux decrement values (Da) of Steidel & Sargent (1987), 
A = 0.0028 ± 0.0004, 7 = 2.46 (fixed). The errors pro- 
vided correspond to the equivalent Icr Gaussian errors, as 
described above. 

The results arc tabulated in Table Bl. All the estimates 
are found to be in accord, including those of Press ct al. This 
is noteworthy in that the estimates that had been made us- 
ing the best value for fa from Press et al. generally gave val- 
ues significantly higher than measured using high resolution 
spectra. Estimating the mean flux including the errors as 
we have done shows that these results are actually in agree- 
ment with the estimates from the higher resolution data. It 
is also noteworthy that the estimates based on the moder- 
ate resolution data of Steidel & Sargent are nearly identical 
to those using the fa fit of Kim ct al. , which is based on 
a compilation of several high resolution spectra. The agree- 
ment in the determinations of the mean flux between the 
moderate and high resolution spectra suggests that moder- 
ate resolution spectra, like those of the Sloan Digital Sky 
Survey, should, in principle, be adequate for measuring the 
mean flux. At attempt to measure the mean flux from the 
Sloan QSO survey was made by Bernardi et al. (2003). The 
values found ai z = 2.41, 2.75, 3.0 and 3.89 are, respec- 
tively (as given in their Fig. 21), (exp(— r)) = 0.74 ± 0.03, 
0.71 ± 0.01, 0.65 ± 0.01, and 0.36 ± 0.02. These all lie sys- 
tematically low compared with the values listed in Table Bl. 
The reasons for the discrepancy are unclear. 

Since the mean flux is the quantity of interest for nor- 
malising the Lya forest simulations, wc use direct measure- 
ments of the mean flux to normalise our simulations rather 
than estimates of the mean Lya optical depth. Historically, 
fa was of interest in the context of models for which the ab- 
sorption was produced by discrete aborption clouds rather 
than modulations in a continuous medium. For the purposes 
of studies of the IGM under the current continuum model of 
its structure, fa has lost most of its significance; the mean 
flux is almost always to be preferred. An exception is in es- 
timating the contribution of known discrete systems to the 
mean flux. Of particular interest is the contribution from 
Damped Lya Absorbers (DLAs). The contribution of dis- 
crete absorption systems to fa is given by (eg, Meiksin & 
Madau 1993; PRS) 



where wi^Nui) is the equivalent width corresponding to a 
system with H i column density A^hi, and At, is the (rest) 
wavelength of the Lya transition. Here, d'^N/{dN^idz) is 
the H I column density distribution at redshift z, assumed 
to vary slowly over the width of an absorption feature. The 
A'hi distribution is estimated to be a power law cx N^^ with 
(3 = 1.5, with a possible high A'hi cut-ofF (Storrie-Lombardi, 
Irwin &: McMahon 1996). Since DLAs are on the square- root 
part of the curve of growth (w oc N^-^), tdla is logarith- 
mically divergent at the upper end. Even if the distribution 
steepens at the upper end, as suggested by Storrie-Lombardi 
et al. , the statistics are so poor that the contribution of 
DLAs to fa is highly uncertain. Seljak ct al. claim to have 
adjusted their error estimate of fa to allow for DLAs, al- 
though they give no details of their procedure. In view of 
the uncertainty in tdla, the most reliable approach is to 
use only spectral regions that are clean of DLA contaminar 
tion to estimate the mean flux. 
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